%% Analyzing the outcome of SEAWAT; genereating graphical output
%
close all; clear all; clc;

load basename
load underneath

movie=0;
%% preperation, read datafiles
if 1
    C   = readMT3D('MT3D001.UCN');
    C   = maskHC(C,[0 1],[0 1]);
    for i   = 1:length(C)   % [kg/m3] convert concentrationvalues to added-salt-content
        C(i).values     = C(i).values.*29;  % as stated in the spreadsheet, saltwater adds 29 [kg/m3]; if changed sere clim
    end
end

if 1
    B_  = readBud([basename,'.BGT']);
    B=mf_Psi(B_);
    zonebudget(B)
end

if 1
    H   = readDat([basename,'.hds']);
end

%% plot well concentration
for i = 1:length(well)     %not tested for multiple well structures
    well    = setWellCout(well,C,i);
    f = figure();
    hold on
    subplot(2,1,1); well.plotCout(i);   title('amount of salt [kg/m3]')
    subplot(2,1,2); well.plotQ;         title('Quantities [m3/day]')
    drawnow;
end

%% plot grid
fig = figure('visible','on');
hold on;

gr.plotMesh;
view([0 -1 0])
well.plot3D('LineWidth',20);

%  plot properties
xlabel('x [m]');                 % set x-label
ylabel('y [m]');                 % set y-label
zlabel('z [m]');                 % set z-label
c=colorbar('horiz');
ylabel(c,'salt [kg/m3]')
axis vis3d equal tight;

%% Movie
if movie
    aviobj=VideoWriter('Upconing');   % create different video for each serie
    aviobj.FrameRate = 1;
    aviobj.Quality = 100; %1=Lowest quality, 100=Max quality
    open(aviobj);
end

%% Plot data
for i = 1:length(C) % for each stressperiod
%     for i=[1 2 length(C)]
    if 1
        if i == 1
            isorange=unique([linspace(-50,50,50),0]);
            
            U.X     = gr.XM;
            U.Y     = gr.YM;
            U.Z     = gr.ZM;
            
        else
            for ii=1:length(isorange), delete(hpatch(ii)); end
            delete(findobj(gca,'Type','line'));
        end
        
        V.X     = flipdim(cumsum(B(i).term{2},3),3);
        V.Y     = cumsum(B(i).term{3},3);
        V.Z     = cumsum(-B(i).term{4},3);
        
        %	frontal view, over the y ax; will display flow from left to right
        for ii=1:length(isorange)
            [f,v] = isosurface(U.X,U.Y,U.Z,V.X,isorange(ii));
            hpatch(ii) = patch('Faces',f,'Vertices',v,'FaceAlpha',0.35,'FaceColor','interp','EdgeColor','none');
            isocolors(U.X,U.Y,U.Z,C(i).values,hpatch(ii));
            set(gca,'cLim',[0 29])
        end
        
        % To plot concentrationslices which could serve as an backgroud for
        % the streamline above
        %     if  i == 1
        %         [sliceH]    = slice(gr.XM, gr.YM, gr.ZM,C(i).values,[],[gr.ym(end)],[],'cubic');
        %         set(sliceH,'FaceColor','interp','FaceLighting','gouraud','EdgeColor','none','DiffuseStrength',1)
        %     else
        %         cdata = sliceCData(gr.XM, gr.YM, gr.ZM,C(i).values,[],[gr.ym(1)],[],'cubic');
        %         set(sliceH,'CData',cdata)
        % 
        %     end
        %     set(gca,'cLim',[0 29])
    end
    
    fprintf('time plotted is %f\n',perdata.time(i));
    
    set(gca,'DrawMode','fast','cLim',[0 29])
    if movie
        h=getframe;
        writeVideo(aviobj,h);
    else
        drawnow;
    end
end
